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We present a formalism for estimating the angular power spectrum of extragalactic gamma-rays 
produced by dark matter annihilating with any general velocity-dependent cross section. The rele- 
vant density and velocity distribution of dark matter is modeled as an ensemble of smooth, universal, 
rigid, disjoint, spherical halos with distribution and universal properties constrained by simulation 
data. We apply this formalism to theories of dark matter with p-wave annihilation, for which the 
relative- velocity- weighted annihilation cross section is av = a + bv 2 . We determine that this signifi- 
cantly increases the gamma-ray power if b/a> 10 6 . The effect of p-wave annihilation on the angular 
power spectrum is very similar for the sample of particle physics models we explored, suggesting 
that the important effect for a given b/a is largely determined by the cosmic dark matter distri- 
bution. If the dark matter relic from strong p-wave theories is thermally produced, the intensities 
of annihilation gamma-rays are strongly p-wave suppressed, making them difficult to observe. If 
an angular power spectrum consistent with a strong p-wave were to be observed, it would likely 
indicate non-thermal production of dark matter in the early Universe. 
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I. INTRODUCTION 

Based on one simple interpretation of astrophysical observations, in the context of ACDM cosmology, it is estimated 
that about 83% of the matter in the Universe is dark matter, and that this matter accounts for 23% of the Universe's 
total energy content 0, [|J ■ One theory that accounts for the presence of this matter is that of the weakly interacting 
massive particle (WIMP). In this paradigm, the WIMP is a new stable particle that is produced spontaneously in 
the early Universe during the Big Bang. WIMP interactions with the Big Bang plasma, for example through WIMP 
pair production and annihilation, keep its abundance in thermal equilibrium until the Universe becomes too cool 
to produce new WIMP particles. Annihilation of these particles becomes rare once the rate of expansion of the 
Universe exceeds the rate of particle annihilation, and the remaining WIMP abundance is said to freeze out. This 
thermal production of a dark matter relic generates the correct amount of dark matter in our Universe if the WIMP's 
rclativc-vclocity-wcightcd annihilation cross section is of the average magnitude \av]f ~ 3 x 10~ 26 cm 3 /s at the time 
of freeze out. If this is the correct theory of dark matter, then we would expect annihilations to be occurring today, 
predominantly in the densest regions of the Universe. Observation of products from these annihilations not only 
would give us information about the particle physics nature of the WIMP, but properties of an extragalactic signal 
also would be rich in information about the large scale structure of matter. 

There is an ongoing endeavor to search for signatures of dark matter annihilation in cosmic signals including gamma 
rays, cosmic rays, and neutrinos. These are looked for: in nearby point sources like the sun, galactic center, and nearby 
dwarf galaxies; in the diffuse galactic halo: and in the extragalactic distribution Q. Indirect signals have already 
indicated unexpected features. PAMELA [4} observes a larger than expected positron fraction in the energy range of 
60 — 100 GeV, and FGST sees more cosmic electrons than expected at around 500 GeV p. It is possible that these 
anomalies will be understood in terms of improved models of emission from supernova remnants [6J] , or pulsar wind 
nebulae Q • Using observations from one indirect signal to constrain these astrophysical models generates predictions 
for other indirect signals [8]. As our understanding of these more standard astrophysical emission processes improves, 
it becomes more likely that emissions from dark matter annihilation might be extracted. If such a signal is to be 
identified, precise theoretical predictions of its properties are imperative. 

Early estimates of gamma-ray mean intensity and angular power spectrum from extragalactic dark matter an- 
nihilation used the spherical halo model of large scale structure [9j. The simplest WIMP model was used with 
av ~ 3 x 10 -26 cm 3 /s and a parametrization of the annihilation spectrum motivated from the minimally supersym- 
metric standard model (MSSM). This formalism was recently generalized to take into account any theory of dark 
matter annihilation, and a study of different particle physics effects on the mean intensity spectrum of annihilation 
was presented 10]. It is found that in many models the annihilation cross-section is velocity dependent and this has 
a large impact on the calculation of intensity. Examples of velocity-dependent effects in the annihilation cross section 
include a p-wave component [ll|, Sommerfeld enhancements and resonances [lS^ . Breit-Wigner resonances fli^ . and 
combinations thereof. In this work, we revisit this general formalism, presenting it in a simpler form, and we extend it 
to the application of calculating the angular power spectrum of the extragalactic annihilation gamma-rays for general 
velocity dependence of the annihilation cross-section. The present work applies this formulation to the case of p-wave 
annihilation (the formalism can be applied to the other cases of velocity-dependent annihilation in future work), and 
offers some preliminary results. 

The halo model of large scale structure seems to be an appropriate paradigm for these calculations. Annihilation 
within smooth halos is dominated in the cores of the halos where the number density is largest. Since halos are 
predicted by simulations to contain dense substructures, these will also need to be accounted for in order to produce 
realistic predictions. Current estimates show that the contribution of substructure to extragalactic annihilation within 
a large halo can increase the signal by a factor on the order of 100, while the galactic signal seen from within the halo 
is increased by substructure by a factor of only a few [l4|. This subhalo effect is not accounted for in this early work 
and will require attention. 

For simplicity, this work assumes that dark matter is distributed throughout the universe in spherical halos. Al- 
though halos in general are predicted by simulations to be tri-axial, their cores are nearly spherical. We assume 
universal radial profiles of the halos' matter density and velocity dispersion, dependent only on the halo's mass and 
rcdshift. The velocity distribution is currently approximated to be isotropic (equal radial and transverse velocity 
dispersion), which is indicated by simulations to be correct deep in the halo cores Where necessary, we assume 
a locally Maxwell-Boltzmann distribution of the particles, specified by the velocity dispersion at each position 
This knowledge is used to determine the average local relative velocity between any two dark matter particles at a 
particular position. All other needed halo properties, such as concentration, are uniformly taken to be at the ensemble 
average for the given redshift and halo mass. 

For this calculation, it also makes sense to use the rigid halo approximation: far from the halo centers, the dark 
matter density is low and annihilations are rare, so we may assume the density vanishes beyond some appropriate 
radius from the halo. Contributions due to overlapping (i.e. merging or unrelaxed) halos are expected to be small 
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and so it is assumed no two halos overlap. In fact, simulations show small halo cores remain relatively intact during 
mergers with larger halos and contribute as substructures within the parent halo. With all of these considerations, 
our model of large scale structure is precisely an ensemble of spherical, smooth, disjoint, rigid, mass-universal halos. 

A comment about notation: there are two kinds of averages that appear frequently in this paper. For clarity, 
averages over all space at a fixed time (or over a shell at fixed redshift) are assumed equivalent to a halo ensemble 
average at that redshift and will be denoted with angle brackets (•) (z). This is in contrast to an average over a 
distribution at a single position (or over an infinitesimal volume) which will be denoted with an overbar • (r). 

This paper is arranged as follows. Section |TT] presents the general expressions for calculation of the mean intensity 
and angular power spectrum of extragalactic gamma-rays from dark matter annihilation. To clarify these results, 
details leading to the formulae in this section are provided in appendices. The analysis of p-wave annihilation is 
discussed in Section IIIII and results from calculations presented in Section IIVI We share our conclusions and outlook 
in Section |Vl Appendix [A] reviews the intensity of gamma-rays from annihilating dark matter and introduces our 
notations. The elements of large scale structure that we require are described in Appendix [B] and are applied to 
the mean intensity and angular power spectrum of the annihilation gamma-ray intensity in Appendix [C] Finally, our 
evaluations of the angular power spectrum require efficient numerical evaluation of particular Fourier transforms of 
halo profiles. We share our algorithms for doing so in Appendix [Dl 



II. INTENSITY FROM POSITION-DEPENDENT ANNIHILATION CROSS SECTIONS 



If we look out with a gamma-ray telescope in direction n along a line-of-sight light shell, with distance and time 
of photon emission specified by redshift z, the intensity of gamma-rays of energy Ej that are due to annihilation of 
dark matter is (see Appendix [A) 

/ 7 (£ 7 ,n) = J -^- ) W((l + z)E 1 ,z)[p 2 av]{n,z) (1) 



where we define 

W( Ej ,z) ee - \ 1 ^L(E 7 )e-^ (2) 
8irmf) M (1 + z) A dE 7 

as the intensity window function. In these expressions, H (z) is the Hubble function, p(h, z) is the dark matter density 
at the specified position, and av(h, z) is the locally-averaged, relative- velocity-weighted, dark matter annihilation 
cross section at that position. Appearing in the window function is the dark matter particle mass todm, gamma-ray 
spectrum per annihilation < ^ L {E 1 ) 1 and the cosmic opacity to gamma-rays t(£" 7 , z) [l7| . 

In s-wave dominated dark matter models where av is a constant, the intensity distribution goes with the square 
particle number density. However, if the magnitude of the annihilation rate varies over the velocity distribution of 
cosmic dark matter, then the appropriate position-dependent field that determines the intensity distribution is p 2 av. 

Given statistical properties of large scale structure gleaned from simulations in the context of the spherical halo 
model (described in Appendix |B|). we can derive corresponding statistical properties of the extragalactic annihilation 
gamma-ray intensity. The mean intensity is simply 

(Ly)(EJ = J ^- ) W((l + z)E, n z)(p 2 Wv)(z). (3) 

where the mean field is found through an ensemble average over dark matter halos. 

(p 2 av)(z) = J d 3 KdM-^(M, z) p 2 h (R\M, z) [av] h (R\M, z) (4) 

Here, the required elements of the halo model are the halo mass function ^p-(M, z), the universal halo density profile 
Ph(R\M, z), and a universal halo profile of the weighted annihilation cross section [av]h(R\M, z) (see Appendix [Bl for 
details of these elements) . 

The angular power spectrum of the intensity signal (detailed in Appendix [C]) is 
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where the relevant power spectrum of p 2 av is simply 

P p2 ^(k,z) = J dM-^(M,z)[?7{lp 2 av} h }(k\M,z) 



dM^-(M, z)b(M, z)JT{[p 2 av} h }(k\M, z) 
dM 



flin(M" 



(6) 



Here 



r(z) 



o H(z) 



is the distance to a position of redshift z, b(M,z) is the halo bias function, 3^7{[p 2 av]h} is the Fourier transform of 
the halo profile, and Pn n (k, z) refers to the linear power spectrum at redshift z. 

The expressions in this section may be applied to any general velocity-dependent annihilation cross section. We will 
discuss the case of p-wave annihilation now, which commonly appears in various supersymmetric models, for example. 



III. APPLICATION TO ANNIHILATION WITH P-WAVE 



For s-wave annihilation, av = [avjo, a constant. Then the intensity spectrum is simply 

dz 



(J 7 ) (£ 7 ) = [Ho 



where 



and the angular power spectrum reduces to 



H(z) 



W((l + z)Ry,z) (p 2 )(z) 



d 3 rdM^(M,z)p 2 (r\M,z), 



P(IJ a {EjJ H(z) 



-W 2 {{1 + z)E yi z) k 2 P p 2 tP i(k, z) 



fc=- 



with 



(k,z) 



dM^(M 7 z)[k?7{p 2 }(k\M 7 z)Y 



dM^(M,z) [k?7{p 2 }(k\M,z)] 



Plin(k,z). 



The quantity k3 r 7{p 2 l }(k\M, z) for the NFW halo profile that we use approaches a constant in the asymptotic k — > oo 
limit (see Appendix ID II) . Note that, due to the normalization with mean intensity, the angular power spectrum does 
not depend on the value of the annihilation cross section, [crujo- In fact, it is a desirable property of the angular power 
spectrum that it is independent of any uniform constants appearing in the intensity distribution, including constant 
intensity boost factors that may be associated with halo substructures or non-thermal relic effects, or intensity 
suppression factors due to p-wave suppression or co-annihilations during freeze out. 
For p-wave annihilation, the velocity-weighted annihilation cross section is 

av = a + bv 2 = [av]o ( 1 H — V 
V a 

where [<rv] — a and b are constants, and the cross section halo profile is simply given by Eq. (|B2|) . In this case, if 
there is significant dark matter annihilation with square relative velocities > a/b, then the distribution of produced 
gamma-rays is coupled to the cosmic dark matter velocity distribution in such a way that regions of high-velocity 
particles will appear brighter. The intensity spectrum with p-wave annihilation is 



<J 7 ) (£ 7 ) = Mo 



H(z) 



W{{\ + z)E 1 ,z)(p 2 [! + —<: 



(7) 
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where 



d 3 vdM^(M,z)p 2 h (r\M,z) 
Xb 



l + —vlh(r\M,z) 
a 



(8) 



= (p 2 )(z) + ^(A2)W 



In these expressions, a 2 denotes the velocity variance of the dark matter at each position, and a 2 h (r\M, z) is the 
associated universal halo profile. In this work, we considered local velocity distributions where the mean square 
relative velocity at each position is related to the velocity variance by v 2 = Xa 2 for some constant A. For Maxwell- 
Boltzmann distributions, A = 6. 

The effects of the p-wave on the shape of the annihilation spectrum are encoded in the relative contribution of the 
new second term, due to the p-wave, given by 



with 



{I^jE^av = a + bv 2 ) Xb 
{l.y) (E 1 \<TV = a) a 



I 1 ^W((l + z)E 1 ,z)(p 2 a 2 u )(z) 
J 1 ^W((l + z)E 1 ,z)(p 2 )(z) ■ 



(9) 



(10) 



Other than the dependence on large scale structure in the ensemble averages, A/ depends only on the details of the 
annihilation spectrum and opacity effects. Note the relative change in intensity diverges for vanishing [av]o since the 
s-wave intensity is zero in this limit. 

The angular power spectrum with p-wave annihilations is 



Ho 



dz 



Cl[El) = t 2 ^f{E~) J H(zj W ^ + z ) E ^ k ^(1+^2)^) 
where the power spectrum is 

dn 



k=- 



(11) 



^ p 2( 1+ M CT 2)(M) 



dM— (M, z) 
dM V ' ' 



J dM^(M, z)b{M, z)37 [p 2 + ™p 2 a 2 uh \(k\M, z 



P\in(k,z) 



P 2 2 {k,z) + 2 — P 2 ,p2 CT 2(fc, Z) + 

a 



Pp 2 <?l,p 2 <?l(k, z). 



(12) 



(13) 



For clarification, the mixed power spectrum is 



2ff 2 (k,z) = (p 2 )(z) (p 2 <7 2 ,)(z)P p 2 7p 2 (7 2(k, z) 



= I dM—(M,z)TJ{pi}(k\M,z)TJ{p 2 h a 2 uh }(k\M,z) 



I (All ( (i77 

J dM— (M, z)b(M, z)37{p 2 h }{k\M, z) J dM— (M, z)b(M, z)37{p 2 a 2 uh }(k\M, z) 



P hn (k,z). 



The biggest challenge in evaluating these expressions is the efficient evaluation of the Fourier transforms. Numerical 
integration of the Fourier transforms for each integrand sampling during the halo mass and redshift integrations is 
more time-intensive than is reasonable. See Appendix [D] for the efficient algorithms we implemented for evaluation 
of these transforms for the case of NFW halo profiles. 

The relative contribution of the quadratic term in av to the angular power spectrum is 



C e {E,\av = a + bv 2 ) 1 + ^Ag(£ 7 ) + 



A 2 ) 



C 0> e(EJav = a) 



[1 + ¥A,(25 7 )]' 



(14) 
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FIG. 1: Left: Intensity spectrum of diffuse gamma-rays from extragalactic annihilating dark matter for the particle 
model described in the text and for matter distributed according to the spherical halo model, and assuming a thermal 
dark matter relic. Right: The relative effect of p-wave annihilation on the spectral shape for this model, per Xb/a. 
Its magnitude of 10~ 8 is typical for any particle physics model. 



where each multipole £ has its own set of power spectrum coefficients 

A c (^7) = = 



k=t/r(z) 



k=£/r(z) 



k=i/r(z) 



It is more convenient to re-express the p-wave effect as 

Ci{E y \av = a + bv 2 ) 



k=t/r(z) 



Co,i{E 7 \av = a) 



= 1 



' Ab' 



where 



Ag(15 7 ) = Ag>(iS 7 )-2A,(i5 7 ), 



A^(E y )-A 2 (E 7 ). 



(15) 



(16) 



(17) 

(18) 
(19) 



It is interesting to note that this has a well-defined finite limit in the vanishing a limit, and that Ajv* does not 
contribute in that limit. 



IV. SAMPLE CALCULATIONS 



We calculated the angular power spectrum for particle physics models with p-wave annihilation components that 
were used in jlOj. We find the results to be nearly universal against the different particle physics phenomenologies, 
suggesting that the effect of p-wave annihilation on the angular power spectrum is largely determined by details of 
large scale structure. We will present results of the calculation for the most typical model. The angular power spectra 
results for the other models were within 50%. 

The typical particle physics model is taken from minimal supergravity (mSUGRA), with parameters m = 
2569 GeV, 7774/2 = 395 GeV, tan/3 = 10, A = 0, and fi > 0. The dark matter is the lightest neutralino Xi 
with mass m^o — 150 GeV. The needed particle physics details of this model were calculated using DarkSUSY 5.0.5 
interfaced with ISAJET 7.78 and FeynHiggs 2.6.5.1 [20|. This model was calculated to have a thermal relic 
density of JIdm^ 2 = 0.114. This chosen example is from the focus point region of parameter space. 

In Figure [TJ we plot the mean intensity spectrum of the diffuse extragalactic annihilation gamma-rays, and the 
relative effect of the p-wave on the spectral shape. If we had chosen a model in the stau neutralino coannihilation 
region of parameter space at low tan /3, where the s-wave of the dark matter annihilation cross section is strongly 
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FIG. 3: The angular power spectrum of extragalactic, diffuse gamma-rays from dark matter annihilation with different 
p-wave components. 



helicity suppressed making the p-wave component strong, we would find that the intensity curve is similarly shaped 
but is as suppressed as the cross-section's s-wave. However, Aj is very similar between the different models. 

To see the effects of the p-wave on the angular power spectrum at the peak energy of the mean intensity spectrum, 
we plot the Ac, coefficients in Figure [21 Here, we have chosen the same model as we mentioned above. If we had 
chosen the stau neutralino coannihlation scenario, we get a similar picture. These coefficients were nearly universal 
for the various particle physics models we explored, varying at most by about 50% from the model shown here. Based 
on the magnitude of the coefficients, we again find that a p-wave will only have a significant effect on the angular 
power spectrum if Xb/a > 10 7 ; that is, if b/a > 10 6 . Unfortunately, the p-wave suppression of these thermal relic 
theories is so large that it makes it unlikely to observe such a model via indirect detection of extragalactic gamma-rays 
[lClj |. It may, however, be possible to consider non-thermal relics with both significant s-wave components and strong 

p-wave strengths. It is interesting to take the general shapes of Ajv 1 and A^ , and put them into Eq. (IT71) for various 
values of Xb/a to see how the angular power spectrum can be affected by the coupling of dark matter annihilation 
to the particle velocity distribution. The results of this exercise are shown in Figure [31 At b = 0, the usual s-wave 
angular power spectrum seen in previous works is reproduced We note that a strong p-wave can significantly 
increase power, more so for large values of I. If a component of gamma-rays of extragalactic origin is determined to 
have an angular power spectrum that is best described by a dark matter annihilation with significant v 2 component 
in its cross section, it would be an interesting challenge to understand the mechanisms that allow such a signal to 
be observable. The magnitude of the effects for the p-wave cross section are very motivating for the consideration of 
other interesting scenarios, such as annihilation resonances at low particle velocities. 
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V. CONCLUSIONS 

Using the formalism of the universal spherical halo model of cosmic dark matter, we have developed techniques 
for estimating the mean intensity and angular power spectrum of gamma-rays due to extragalactic annihilating 
dark matter for general velocity-dependent annihilation cross section. The formalism for these calculations has been 
simplified in the case of the mean intensity of the signal, and is new for the angular power spectrum. It is found that 
the important spatially-dependent field that determines the distribution of intensity is the quantity p 2 av, where p 
is the dark matter density, and 7tv is the locally averaged, relative- velocity-weighted, dark matter annihilation cross 
section. While the mean intensity can be numerically calculated in a straight-forward way, the angular power spectrum 
is numerically challenging. It requires evaluation of the Fourier transform of the universal halo profile of p 2 av, which 
is dependent on the halo's mass M and redshift z. Efficient evaluation of the Fourier transform is necessary in order 
to complete numerical integrations over M and z. This is prohibitively time-consuming if a general-purpose integrator 
is used. We succeeded in developing a sufficient algorithm for this task for the case of rigid NFW halo profiles and 
annihilation cross sections with a p-wave component. 

Thus, we apply our formalism to the specific case of dark matter annihilation with a p-wave component. We 
find the p-wave only affected the angular power spectrum of extragalactic annihilation gamma-rays for big p-wave 
strengths, b/a, of about 10 6 or higher, for our model of large scale structure. Since p-wave suppression of intensities 
of gamma-rays from annihilating thermal relics makes them very difficult to observe, the observation of a component 
of gamma-rays with an angular power spectrum consistent with strong p-wave annihilations would be an indication 
of interesting new early-universe physics, and/or new enhancement mechanisms of the annihilation intensity. 

Having established that velocity-coupled annihilation of cosmic dark matter can have significant effects on the 
angular power spectrum of produced gamma-rays, now it will be interesting to understand the angular power spec- 
tra associated with other realistic features of annihilation cross-sections, especially those which would enhance the 
annihilation intensity. Breit-Wigner resonances near the dark matter particle rest mass would produce interesting 
velocity-dependent effects in the present distribution of cosmic dark matter. Also, Sommcrfcld enhancement of anni- 
hilation, and Sommerfeld resonances result in new shapes of av that favor slow moving particles, as opposed to the 
p-wave, which favors the annihilation in regions where particle are moving rapidly with respect to one another. Our 
general formalism can be applied to all these cases. 

It is likely that halo substructure has a significant effect on the annihilation signal, perhaps increasing the extra- 
galactic mean intensity by a factor as high as 100, as opposed to the intensity from annihilations within our galactic 
halo being increased by a factor of a few. If, for different halos, this intensity factor varies little with halo mass, it 
would not affect the angular power spectrum, but significant halo mass dependence would have interesting effects. It 
is important to add realistic descriptions of the substructure to this formalism to improve its predictive power. There- 
fore, this will be considered in future work. It is also important for further research to understand the robustness of 
these calculations against the uncertainties of the large scale structure, in order to have the ability to extract particle 
information from an observable signal, but as well, to understand how effectively a signal will constrain the properties 
of the large scale structure distribution. 
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Appendix A: The intensity of annihilation gamma-rays 

In a gas of particles (with number density n) that may annihilate with one another, the cross section of annihilation 
a is defined as the rate T p of annihilations per target particle divided by the incident flux on the target nv where v 
is the relative velocity of the incident particle and target particle. The mean annihilation rate per target at a given 
position is T p = nav. The annihilation rate per unit volume at a given position is T = \nT p = ■^n 2 'av. 

The rate of energy emitted due to annihilations in a given density of dark matter is given by the power emissivity: 
the emmitted energy of photons with energy between E 1 and E 1 + dE 1 per unit volume, time, and energy range dEj. 

P 7 (£ 7 )d£ 7 = ^n 2 avE 7 dN 7 

where d7V 7 is the number of number of photons per annihilation with energy between E 1 and E 1 + dE 1 . We call j^f- 
the differential photon spectrum per annihilation. 
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In a flat Friedmann-Robertson- Walker cosmology, using coordinates in the cosmological rest frame, the proper 
volume of space with solid angle dft and thickness dz at redshift z is 



dV = [a{z)dr][a 2 {z)r 2 drt] = r 2 drdSl 

where a is the cosmological scale parameter. 

The luminosity dL 1 {E 1 ) of photons of energy E 1 emitted from this region of space because of annihilations is 

dL 7 (E 7 ,z) = P^(E^z)dV = ±n 2 avE, ( ^(E 7 ) * r 2 drdSl. 

Assuming isotropic emission, the photons emitted by this volume pass with uniform flux density through any sphere 
centered on the source. The sphere on which we sit, centered on the source, has proper surface area 

A = 47rrV(0) = 47rr 2 . 

The total luminosity on this shell (energy of photons emitted from the source with energies between E 1 and E 7 + dE 7 , 
per dE-y, per unit time of emission) is redshifted: the cosmological redshift of photon energy due to the expansion of 
the universe is cancelled by the redshift of the energy bin d-£? 7 ; the arrival rate of photons is redshifted giving one 
factor of (1 + z)^ 1 . Observation of photons of energy E 1 means photons of energy (1 + z)E 1 were emitted. Hence, 
the luminosity of photons on the observer's spherical shell with energy E 7 from the source at redshift z is 

dL^E^z) = dL ^ 1 + Z ^ Z \ -r{(l+z)E y , Z ) 

1 ~\~ z 

where t(E 7 , z) is the opacity of the universe to gamma rays |17| . The photon flux on the sphere, or surface brightness, 
due to a source at position r and redshift z is 

dL'(E~,v,z) 1 1 „ diV-v in i \n \ 

dS 7 (£ 7 ,r,z) = ^ = -^—^n 2 {v)av{v){l + z)E^{{l + z)R,)e^ 1+ * E ^drtil. 

The net specific intensity (number of photons of energy Ej observed per bin di? 7 , per unit time, per source solid 
angle, per normal photon collecting area) is found from a line-of-sight integration in direction n: 

7 7 (£ 7 ,ri) = J = J ■^- ) W((l + z)E„z)[p 2 av}(h,z) (Al) 

where the coordinate distance is related to redshift via dr — —jf^pj with the usual Hubble function given by H(z) = 

Hoy^fl m (l + z) 3 + tt\ in terms of the Hubble constant Hq, local matter abundance fl m , and local dark energy 
abundance Ha = 1 — Q m - The important spatially dependent field p 2 aiJ, where p is the density of dark matter, is 
weighted by the intensity window function 

lf ' e ""° K > ln-'-)- S" i '""' M <A2) 

with ttidm being the mass of the dark matter particle. 



Appendix B: Modeling large scale structure with the spherical halo model 

Statistical properties of the spatial distribution of p and av at each redshift can be modeled with the disjoint 
spherical halo model. Since p throughout the universe is greatest at the cores of halos that are nearly spherical, it is 
reasonable to assume that the dominant contribution to an annihilation signal is due to dark matter consisting of an 
ensemble of disjoint spherical halos, each at position with some mass M,. 

1. The point distribution of halos 

The point distribution of Nh halos at redshift z is 

N h (z) 

p h (r,M,z)= J2 i (3) (r-R,(#(M-M t (z)). 
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The halo mass function is the ensemble average of the point distribution: 

dn 



normalized such that 



/ 



d ^(M,z) = (p h )(M,z) 



(All/ 

dM— M = (p)(z) = Pc n DM (l + z) 3 



where p c = is the local cosmological critical density, and 51dm is the local cold dark matter content. 
The halo overdensity can be defined formally as 

which has zero ensemble average. This is correlated to the matter overdensity 

(p)( z ) 

This correlation is loosely described by the halo bias function 

' 5 h {v,M,z) 



b(M, z) 



S p (r,z) 



which will be described, for our purposes, more precisely below in terms of the halo power spectrum. 

The matter correlation function over 2 positions r x and r 2 at the same redshift z is £(rx, r 2 , z) — (<L(i"ii z)5 p (r 2 , z)). 
It only depends on the separation r = \r 1 — r 2 | of the 2 points: £(r, z) = (S p (rx, z)S p (r 1 + r, z)). The power spectrum 
is the Fourier transform of the correlation function: P PjP (fc, z) = j d 3 re~ jk ' r £(r, z). The evolution equations for S p 
are non-linear. When linearized, the linear solution results in a power spectrum valid on linear scales called the linear 
power spectrum. 

P PiP (k, z) ~ Pii n (fc, z) for small k. 

In analogy, the full halo correlation function is ^(ri, Mi, r 2 , M 2 , z) = (Sh(ri, Mi, z)8h{v 2 ,M 2 , z)) . We can see that 
this function has a singularity at ri = r 2 by expanding 

Z i AT s (Ph{ri,M 1 ,z) Ph (r 2 ,M 2 ,z)) 
£ h (r 1 ,M 1 ,r 2 ,M 2 ,z) = y \An (M ~s 1 

and splitting the point distribution 2-moment into diagonal and non-diagonal pieces. 

/ N h N h \ 

(p h {r 1 ,M 1 ,z) Ph (r 2 ,M 2 , z)) = / ^£ 5{3)( - ri ~ R *)<W - ^)<5 (3) (r 2 - Rj)S(M 2 - M 3 ) \ 

\ i=l j=l / 



]T J2 <* (3) (n - Ri)<*(Mi - M t )S^(r 2 - R j )S(M 2 - M 3 ) 



\ i jj&i I 

+ fa tf 3 >(n - R i )6(M 1 - Mi)S^(r 2 - Ri)S(M 2 - Mi)\ 
= Cj?\r 1 ,M 1 ,r 2 ,M 2 , z) + <^(n - r 2 )<5(M 1 - M 2 )^(M 2 , z) 

(2) 

This expression makes the singularity explicit. Here, the function C h is introduced as the non-diagonal part of the 

(2) 

halo 2-moment. Note that if one restricts to disjoint halo ensembles, then C h — formally. 

The non-diagonal part of the halo correlation function, known simply as the halo correlation function [2l| , is now 
everywhere finite and is defined as 

t ( at at ^ Cff^n.Mj.ra.Ma,-*) , 
thfri, Mi, r 2 , Mo, z) = — f ; 1. 

[ ' ^(M 1 ,z)^(M 2 ,z) 
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Then 

it Ayf M \ c( m u , , ^(n -r 2 )J(M!-M 2 ) 
^(ri,Mi,r 2 ,M 2 ,z) = &(ri, Mi, r 2 , M 2 , z) + ^— — — . 

The halo power spectrum is defined accordingly as 



P h , h (k,M u M 2 ,z) = 




- ik -^ h (r l ,M 1 ,r l + r,M 2 ,z). 



Its correlation to the matter power spectrum is encoded in the halo bias function. For our purposes, the precise 
meaning of the halo bias is through the relation 

P hth (k, Mi, M 2 , z) = b{M u z)b{M 2 , z)P PtP (k, z). 

Although recent works show it is likely that these bias functions have some dependence on scale k (22[, we did not 
consider such models in this work. It turns out that the halo power spectrum is only relevant for us on linear scales 
where P PiP may be substituted with the linear power spectrum. For less abstract and more physical discussion of 
many of the quantities in this section, we recommend Peebles (1980) [23| where the discussion is in terms of galaxies 
instead of dark matter halos, but is still relevant none-the-less. 

For our calculations, the knowledge we require of the halo distribution is the mass function and halo bias function. 
For this work, we used the Sheth-Tormen halo mass function [24| and the halo bias due to Sheth, Mo, and Tormen 
[25| . For the matter distribution, we used the best- fit cosmological parameters from WMAP5 [l[ and the linear power 
spectrum of Eisenstein and Hu [26[ with neutrino free-streaming and gravitational wave effects neglected. Dark matter 
halos were taken to have a minimum mass cutoff of 1O _6 M0. 



2. Universal Halo Profiles 



In this section, we describe radial profiles Fh{r\M, z) of halos that are seen to be universal in large N-body simu- 
lations of self-gravitating, collisionless particles. These radial profiles in our model are written to depend only on the 
halo's mass M and observed redshift z. It is possible that, in the future, certain variations between halos of same 
mass and redshift might be taken into account by extending the universality to depend on other halo properties (such 
as halo concentration, formation redshift, angular momentum, shape parameters, etc.), and appropriately extending 
the halo mass and bias functions to also depend on those halo variables. 

The first hint that the phase space of relaxed dark matter halos is universally stratified came when Navarro, Frenk, 
and White [l?} observed that the spherically-averaged radial density profiles of the halos in their simulations were 
well-described by what is now known as the NFW profile 

for r < i?vir, 



Ph (r\p a ,r„R vif ) = \*:( 1 +*:) ' ' (Bl) 

[0 otherwise. 

Halo models typically truncate the halo profile artificially at a virial radius R V11 (M, z) of the mass M halo located at 
redshift z, which is fine for our calculations. Various mass-radius relations have been used in the literature and the 
effect of many of them on the halo mass function were studied in |28[. Based on the results of that study, we adopted 
the convention 

M = -TTB% il (M,z)A viT {p)(z) 

with A v ; r = 180. Although there are more modern density profiles that better describe the most recent simulations 
[l5j , the NFW profile still provides a reasonable description of simulated halos and is a good place to begin development 
of techniques for numerical calculation because of its simple analytic form. The usual method for writing (|B1|) in its 
universal form is to use our mass-radius relation to define a halo concentration 

i?vir(M,z) 



c(M, z) 



r s (M,z) 



which we take to be distributed according to the model in [29|]. This determines r s (M, z). Solving for M by integrating 
over the halo density requires 

/ \ A vir p(z)c 3 

Ps(C, z) = — r — =r 

3 ln(l +c)-j?- 
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which fixes p s {M,z). Thus, we determine the NFW density profile distribution ph(r\M, z). 

It is supposed that the distribution of particle velocities u at each position r of a relaxed universal halo is well 
described by a velocity distribution function /^, iU (u|r, M, z). The halo's velocity dispersion a u h, observed in simulations 
to be nearly isotropic at the halo cores and taken as such in this work, is fixed by the observation of stratification of 
the so-called pseudo-phase-space-density in simulations according to [30jj 

Ph - a 

—5- oc r . 

Coupled with the radial Jeans equation, which describes the phase space evolution of self-gravitating collisionless 



systems, the Dehnen-McLaughlin (D-M) halo profiles and velocity dispersions are determined [3l|. Their critical, 
NFW- like solution has a = y|, consistent with values of the radial dispersion seen in the simulations 15j. The scale 
radius tq for these profiles is at the position where — dln/?h(ro)/dlnr = 6 — 2a. Matching the NFW profile to the 
critical D-M profile at ro, we find an NFW velocity variance halo profile given by 



where /3=|(a — 1) = |y, the scale variance is 



4/3 



a 2 = 12^^-1 U-\ 



0-1/3 



{l-Pf-'psTl 



and k = for the critical D-M solution. 

What we are specifically interested in for the annihilation cross section is the distribution of particle relative 
velocities v at a given position. If one knows //,.. u (u|r, M, z) then the 2-particle distribution of, say, square relative 
velocities fh ,v 2 {v 2 \r, M, z) is easily determined. In the case where // ljU is Maxwell-Boltzmann, one finds the mean 
square relative velocity is simply 

v 2 = 6<rl, 

a relation that would hold at each position of the halo. Since the velocity distribution is well-established to deviate 
from Maxwell-Boltzmann [32[ , one would expect a generalized relation such as 

^(r\M, z) = X(r\M, z)a 2 uh (r\M, z) 

to hold in relaxed spherical halos. For this work, we will continue by approximating A to be constant and treat it as a 
free parameter of our model, to eventually be determined from simulation data. We used A = 6 for calculations that 
required a specific value. 

Taking the velocity- weighted dark matter annihilation cross section [av)(v 2 ) as a function of v 2 , the halo's mean 
cross section profile is determined 



[av] h {r\M,z) = J d[v z ] [av]{v z )f Kv 2{v z \r,M,z). 

For p-wave annihilation with av = a + bv 2 , this process is trivial. 

[av] h (r\M, z)=a + b\a 2 uh {r\M, z) (B2) 

Appendix C: The mean intensity and angular power spectrum of extragalactic gamma-rays: application of 

large scale structure 

From (|A1|) , the mean intensity of annihilation gamma-rays is found simply from averaging over ensembles of dark 
matter halos 

(J 7 ) (£ 7 ) = J ^)W((1 + z)E„ z) (p 2 av)(z). (CI) 
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In the disjoint halo model, an ensemble of halos at redshift z has 



N h (z) 

[p 2 av}(r,z) = (%(T-'Ri(z)\M i (z),z) [av] h (r - R(z) | M 4 (z), z) 
»=i 

where r are a global set of coordinates at the time associated with redshift z. For disjoint ensembles, at most one 
term contributes to the sum at any given position r. This allows us to express the ensemble average in terms of the 
halo mass function. 

N h (z) 



(p 2 av){z) = I^J d 3 RdMp^(r — H\M,z) [av] h (r - R | M, z) ^ <5 (3) (R - R»(-z)) 5{M - Mi(z))j 

= J d 3 RdM^(M,z) p 2 h (R\M,z) [Wv] h (R\M,z) (C2) 

We explore the angular anisotropics in the intensity signal by determining its angular power spectrum defined as 

C l = {\a tm \ 2 ) 

with spherical harmonic coefficients obtained from 

5j (n,f; 7 ) = ^Hl|z2-l = f; J2 «t m (E 7 )Y tm (n), 



or 



a tm (Ej = <j> dI2Jj(n,£; 7 )y^(n). 

J^jf^ J ^{[p 2 av](n,z)-(p 2 av) (z)}W((l + z)E 7 ,z)Y e * m (n) 



(I 7 )(E. 



where, as usual, 



x _ - p2(JV _ i 

Op 2 (TV / 9 \ -L- 

\p A av) 



Then 



where 



1 /" d7 Hz' 



F t {z,z') = J dndn / (S p 2^(h,z)5 p 2^(h\z'))Y; m (n)Y em (h'). 
We'll see shortly why Fi is independent of to. To simplify it, write in terms of the power spectrum of the p 2 av field, 



where 

rx dz 



J 

Jo 



H(z) 

is the distance to redshift z, and similarly for r'. Applying Rayleigh's formula 

oo i' 
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and the orthogonality of spherical harmonics, one finds 

2 f°° 

F e (z, z') = - dk^P^k, z, z')j e (kr)j e (kr'). 

I" JO 

One would not expect any significant correlation between regions of different rcdshift along a line-of-sight. One way 
this is realized is when P p 2 WV is a slowly-varying function of k. In this case, it is a good approximation to treat it as 
a constant at wave number where je(kr) is maximized. Since je(x) has its maximum near x — £, we can approximate 
the power spectrum by its value at k = £/r(z). Then orthogonality of the spherical Bessel functions 

POO 

^ dkk 2 Je (kr) Je (kr') = —S(r - r') = -^H(z)S{z - z') 



gives 



F t (z, z') « 6(z - z')^L k 2 P p ^(k, z) 



We then finally express the angular power spectrum as 



Q(i?7) " P (l!f(E,) I W) W2{{1 + z )E„z)k*P p ^(k, *) . , , (C3) 



where we denote 

P p 2wv(k,z) = (p 2 av) 2 (z)P p 2 W! j{k 1 z). 

To derive the expression for the power spectrum of p 2 av, consider the correlation function at two points ri, r 2 at the 
same redshift z. 

(S p 2^{r 1 ,z)6 p 2^(r 2 ,z)) = 1 

(p 2 av) (z) 



Recalling 



(p h (Ri, M u z)p h {R 2 , M 2 ,z)) = — (Mi, «) — (M 2 , *)[&(Ri, M u R 2 , M 2 ,z) + 1], 



the 2-moment becomes 

'ym]{T 1 ,z)\p t m}]{T 2 ,z)) = (j2J2^ 2mJ ^ ri -Ri\Mi,z)[p 2 av] h {r 2 -R 3 \M 3 ,z^ 

= j d 3 R 1 dM 1 d 3 R 2 dM 2 [p 2 OTj] /l (r 1 -•R.i\M 1 ,z)\p 2 mj\ h {T 2 -R 2 \M 2 ,z)(j3 h (R u M 1 ,z)p h (R 2 ,M 2 ,z)) 

d 3 R 1 dM 1 d 3 R 2 dM 2 — (Mi,z) — (M 2 ,z)[p 2 av] h (r 1 - R^, z)[p 2 mj] h (r 2 - R 2 \M 2 , z) 

x ^ h (R 1 ,M 1 ,R 2 ,M 2 ,z) + j d 5 RdM^{M,z)[p 2 av] h {r 1 - R\M, z)[p 2 Wv] h (r 2 - R\M , z) + (p 2 Wv) 2 (z) 

where the second term in the last line is due to the singularity in We therefore find the correlation function to be 

dn [p 2 av] h (ri — R\M, z)[p 2 av]h(r 2 — R\M, z) 



{5 p 2^{v u z)5 p 2^{v 2 ,z)) = j d 3 RdM-^(M,z)J 



(p 2 av) (z) 



+ /d'RtdMxd'RadMA 

J dM dM (p 2 av) (z) 

x ^ h (R 1 ,M 1 ,R 2 ,M 2 ,z). 

This simplifies significantly in momentum space. If we determine the Fourier transform of the halo profile 



TJ{[p 2 av] h }{k\M,z) = J d 3 re- 4k ' r [p 2 av] h (r\M, z), 
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the power spectrum can be written 

P p 2^(fc,z) = (p 2 av) 2 {z) / d 3 r(<5 p 2 ?fTJ (r 1 ,z)^2 ?fIJ (r 1 +r,z))e-' ir ' k 



dM^(M, z) J7{[p 2 av} h } 2 (k\M, z) 

(i n c\ u 

dMl dM 2 — (M 1 ,z) — (M 2 , z) ■J%\p 2 av\ h }{k\M ll z) ?7{{p 2 av] h }(k\M 2 , z)P h (k, Mi, M 2) z). 

The first term, the one-halo term, dominates at small scales (large k) and the second term, the two-halo term, 
dominates at the large scales, in the linear regime. So in this expression, it is correct to use 

P h (k,Mi,M 2 ,z) = b(Mi,z)b(M 2 ,z)Plin(k,z). 

Therefore, the power spectrum is simply expressed as 



P p2 ^(k,z)= I dM-^-(M 7 z)[H7{[p 2 av} h ]{k\M,z)] 



, W " i I' ! ' N rrnr , II I: W ' I 

dM 

2 

flin(M). (C4) 



J dM^-(M, z)b(M, z)TJ[[p 2 av] h }(k\M, z) 



Appendix D: Numerical evaluation of needed Fourier transforms of rigid NFW profiles 

Although there exist some very good general integrators for Fourier transforms (33j . their use is not very feasible 
in this calculation. The transforms appear in the integrand of the halo mass integration, and that result is then 
integrated over redshift. The number of evaluations required for precise calculation is very large, and takes too long 
to complete when using a general-purpose integrator. Since these functions are over a 3-dimensional space (k\M,z) 
that stretches over a large range of scales, it is also not feasible to fill a data table for interpolation. 

For the rigid NFW profile, a closed form solution is available for yxjp 2 }, which has allowed efficient calculation 
of s-wave angular power spectra in previous works. No such closed form is available for the non-analytic 5T{^ff^ h }. 
Nevertheless, we were successful in developing an efficient numerical algorithm for efficient evaluation of this func- 
tion, described below. One of the challenges for calculations of angular power spectra of extragalactic dark matter 
annihilation products is the development of efficient numerical methods to evaluate 3 r 7{p 2 l \av)h} for a given model's 
halo profiles and annihilation cross section. This calculation would have taken weeks to complete using the quadpack 
general purpose Fourier transform integrator, qawf. With the algorithm described in this section, the results in this 
paper were able to be evaluated within a few days of run time on a desktop computer. 

1. TJ{p 2 h }(k\M,z) 

This Fourier transform can be expressed as 

__, 2wm „ 2 3 f 2 4 + 3c ll + 15c + 6c 2 -[(l + c)fcr s ] 2 Si(kr s c) 
3V{fi}(k) = ^p 2 s rl \ -- + cos(fcr s c) + L \ — '- sm(Ax s c) 1 



3 6(1 + c) 2 ' ' 6kr s (l + c) 3 y ' kr s 

_ (, (kr 3 ) 2 j^_krA{ cos(fcr s ) sin(fcr s )\ (Ci(kr s (l + c)) - Ci(kr s )\\ . . 

V 6 kr s 2 J \-sm(kr s ) cos(kr s )J \Si(kr s (l + c)) -Si(kr s )J j [ ' 

where Si and Ci are the sine integral and cosine integral, respectively, for which efficient numerical methods for 
evaluation already exist [34j |. Evaluating the line-of-sight integrand for the angular power spectrum near z = 
requires the Fourier transform to be evaluated in the /; H> oo regime. One finds that for kr s 3> 1, 



c(l + c)' 



cos(fcr s c) 



0((fcr s )" 



Unfortunately, in the Bullock, et al. model of halo concentrations, the mean halo concentration vanishes at a maximum 
halo mass scale. To stay true to the definition of the model, we want to be able to evaluate the transform in the 
vanishing concentration regime. Here, one should use 
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If c <C 1 and c <C kR vlr (equivalently, kr s ;§> 1), then 

k 37{p 2 h }(k) = 4Tr(r sPs ) 2 jsi(fci? vir ) - 2[l - cos(fci? vir ) 

+0(c 3 ) + 



kR v 1 I 



+ 3 



Sin(fci? v i r ) — fci? v j r COs(fci? v ir) 



\kR v 



In the case where c< 1 and c ^ kRyiv, then we must have kR v „ <C 1 and can use 

,fcit 



fc yT{^}(fc) = 47r(r sPs ) 



2 a- j tvir 



1 / 1 

3 V ~ (1 + c) 3 J 18(1 + c) 



(fci? vir ) 2 + ((fci? vir ) z 

1 + c) V 



2. 3T{p^J(fc|M ( 2) 
This Fourier transform is simply expressed in the form 

L-n-T- 2 n 2 + 

-S(kR viv ,c) 



TT{pW uh }(k) = 



where we have defined 



S(x, c) 



^(1+*)* 



(D2) 



(D3) 



with (3 = 17/27 as previously defined, and q = 16/3 for the NFW profile. The important result that allows efficient 
evaluation of §(a;, c) for a wide range of scales for x and c is the set of expansions (see Appendix ID 3p 



§(x, c) 



c H \ - (g)p 
(/3) 3 [C/ (/?,/? 



•g+l,-i- 



iJlGSjjS+p+ljia;) 



1 + c 



c 



(l + c)« ^ 



2(?) p 9f[e te tr( P + 1,^-? + l,-ix) 



1 + c 



c < ct 



p 

, C > C T 



(1)1) 

where ct is an appropriate transition concentration. The truncation errors of the two expressions were found to be 
of the same magnitude near c = 0.8, making it a reasonable value for ct- Also in the expression appears the gamma 
function F(x), the Pochhammer symbol 

(q) p = ^+2)-= q ( q + l)( q + 2)...( q +p-l), 

the confluent hypergeometric function of the first kind i-Fi(a; b; z) (expressed in the notation of a generalized hyper- 
geometric function), and the confluent hypergeometric function of the second kind U(a, b, z). 

For c < ct, if x is small (say < 4), then the hypergeometric functions are most efficiently evaluated with their 
power series 



1 F 1 ( i 9; J 9+p+ 1;kb)1 = £ 



2n+l 



(-l)"(/3) 2n+1 x 



For larger values of x, the functions are quickly determined from the recurrence relation 

(3 + P ' 



or 



1 F 1 ( ; 3; ; 3+p + l;iar) 

iFi(/9;/9+p+l;ia!) 
iFi(/9;/9+p+l;ia;) 



P 



I -i P + P - ) 1F1 (/?;/?+ p; is) + i - + P 1 1F1 (/?; (3 + p - 1; ix) 



_ p+ 




3? 


p 






_p + 






p 







X 

iFi(/3; P + p; ix) 
iFi(/3;/3+p; ix) 



/3 + P+l ( 

a; 

/3 + p+ 1 



iFi(/9; /3 +p; is) - iFi(# f3+p~l;ix) 
iFi(/9; /3 +p; is) - iFi(# £ +p - l;iar) 
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Since i-Fi(/3; f3; ix) — e lx , then we only need the numerical evaluation of i-Fi(/3; /3 + 1; ix) to be able to determine the 
rest of the sum's hypergeometric functions using the recurrence relation. The power series is suitable for x < 10: 

lFl (P;P + l;ix) = f^-^-^. 

n— 

We find the asymptotic expansion converges appropriately for x > 27: 

1 F 1 ((3;(3 + l;ix)~T(!3 + l)exp(i^)x-e + /3 ^ (n - /?)„_i exp 



n=l,2,3,. 



(*-t) 



For 10 < x < 27, these series' do not converge sufficiently with double machine precision arithmetic. For this short 
range of x, it is not too much of a burden to evaluate the function via numerical integration 



1 F 1 (p;p + l;ix)=p [ e^t^dt. 
Jo 



For large concentrations c > ct, there are two components. The first term depends only on the ratio x = x/c and 
requires the evaluation of $s[U((3, (3 — q + l, — ix)]. We use the perturbative expansion for x < 5, for which a convenient 
expression is 



?[t/(/3, p-q + 1, -ix)] = £ \{-l)W\n mod 2)^-^-^^ 



n=0 



+ ; 



T(q) (p-q + l) n 
(-1) L " /2J (q)n 



2T((3)T(q ~ P + 1) C§„ (9 - £ + l)n ' 

where we have introduced: the modulo 2 operation 



,9-/3 



n! 



0, n even, 

1, n odd, 



n mod 2 

the floor operation [rcj being the largest integer < x, and we define the trigonometric function 

GS n (x) 

The asymptotic expansion 



cosx, n even, 
sin a;, n odd. 



n=0,l,2,... 



n! 



works sufficiently for x > 40. For the little remaining range of 5 < a; < 40, we simply numerically evaluate the integral 
representation 



r(/3)Z[U(p,p-q + l,-ix)] 
To evaluate the functions in the sum, 



sin(xt) 



dt. 



e ix f/(p + l,/3-g + l,-ia;) 



= sinx 5ft 



C/(p+l,/3-g + l,-ix) 



+ cos X 9 



LT(p + 1, /3 - q + 1, -iar) 



we again can make use of recursion relations 



5ft 



U(p+l,/3-q+l,-ix) 



)] ^{x3[c/(p,/3 

J p(p + q-p) L L 



g+l.-ix) + (2p + <? - ,3 - l)5ft LT(p,0-(j+l,-iaO 



-3? 



tf(p- 1,0 -«+!, -»*)]} 



(p + 1, /3 - g + 1, -ix)] = -— -L— ^ {-x$t[u(p, p-q+1, -ix)] + (2p + q-p-l)Z [u(p, 0-q + l, -ix) 



U(p-l,f3-q+l,-ix)]}. 
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Since f7(0, • , • ) = 1, we only require the evaluation of U(l, /3 — q + 1, — ix). For x < 4, 



U(l,/3-q+l,-ix) = 



q-/3 



E 



n=0 



(_1)«/2(1 - „ mod 2) Trx 9- " (_i)L(«+i)/2j 



— y 



n=0 



(/3-g+l)„ 2T(?-fln!eS B+1 (^) 

(_l)(n+l)/2(„ mod 2 ) tto;?-^ (_1)L«/2J 



03 -?+!)„ 2r(g-/3) n!e§n (^^ 



was used, and for x > 45, we evaluated 

[/•(l,/3-g+l,-ix) 



£ (g-/3 + i)„^)-(" +1 ). 

n=0,l,2,... 



For the mid- values of x, we numerically integrated 

17(1,0-3 + 1, -is) = 



dt. 



For very large values of x, the recursion relations will fail, due to loss of precision from subtracted quantities being 
very near each other. In this regime, we can evaluate the hypergeometric function in each term of the sum from the 
asymptotic series 



(? x U{p + l,p-q+l,-ix)\ ~ J2 (-l) l{n ~ p)/2l eS n+I (x)(p+l) n (q-f3+p+l) r 



-(n+p+l) 



n=0.1,2, 



3. Derivation of Equation (|D4[) 



We begin with the case of c < ct — 0.8 by expanding (1 + 1) q in Equation (|D3|) as a power series, and rescaling 
t — > xt/c to get 



m=0 



where 



I n (x) 



t n siat dt. 



Then we let k = c/(1 + c), write the expression in the form 



S(x, c) 



(l + c)i ^ [ ' T(q)m\ xP+ m (1-k 

m— 



X«) 



(1 - n)i+ m ' 



and expand the k expression in a power series with shifted indices 

fi - (<z+m) - V r (g + w + p) r:P = v 

^ T(g + m)p! ^ 



r(q + P ) 



p=0 



Swap the order of summation to find 



b(x,C) = — r— > — — -r 

+ r(g) 



E 



— ' T(q + m)(p — m)\ 

p—rn w ' v± 7 



(-i) m ^ +m -i(x) 



i\(jp — m)\ xP +r 



1 + c 



after substituting c back into k. The quantity in the square brackets can be rewritten using the integral representation 
of the confluent hypergeometric function 
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(convergent for $t(b) > 3?(a) > 0), giving us 
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where, in the second line, we rescaled t — > t/x. 

For the case of large concentrations c > ct , we break § into two terms 

S(x, c) = §i ^— ^ — §2(2;, c) 
sin(it) 



with 
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The first term is simply 



Si(S) = 3 



-dt 
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given that the confluent hypergeometric function of the second kind has the integral representation 
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if 5R(a) > and > 0. For the second term, first substitute t — > t — c 
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and then substitute t — > t/c 
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with k = c/(l + c), as before. As is appropriate for large values of c, we expand as a power series about 
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